{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Changing the input current when solving PyBaMM models\n",
    "\n",
    "This notebook shows you how to change the input current when solving PyBaMM models. It also explains how to load in current data from a file, and how to add a user-defined current function. For more examples of different drive cycles see [here](https://github.com/pybamm-team/pybamm-data/releases/tag/v1.0.0).\n",
    "\n",
    "### Table of Contents\n",
    "1. Constant current\n",
    "1. Loading in current data\n",
    "1. Adding your own current function"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant current  <a name=\"constant\"></a>\n",
    "\n",
    "In this notebook we will use the SPM as the example model, and change the input current from the default option. If you are not familiar with running a model in PyBaMM, please see [this](../models/SPM.ipynb) notebook for more details.\n",
    "\n",
    "In PyBaMM, the current function is set using the parameter \"Current function [A]\". Below we load the SPM with the default parameters, and then change the the current function to be an input parameter, so that we can change it easily later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import pybamm\n",
    "\n",
    "os.chdir(pybamm.__path__[0] + \"/..\")\n",
    "\n",
    "# create the model\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# set the default model parameters\n",
    "param = model.default_parameter_values\n",
    "\n",
    "# change the current function to be an input parameter\n",
    "param[\"Current function [A]\"] = \"[input]\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can set up a simulation in the usual way, making sure we pass in our updated parameters. We choose to solve with a 1.6A current. In order to do this we must pass a dictionary of inputs whose keys are the parameter names and values are the values we want to use for that call to solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "174048e3ccf74930a39a249dc0723975",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=600.0, step=6.0), Output()), _dom_classes=('…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x79ada1f321e0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set up simlation\n",
    "simulation = pybamm.Simulation(model, parameter_values=param)\n",
    "\n",
    "# solve the model at the given time points, passing the current as an input\n",
    "t_eval = np.linspace(0, 600, 300)\n",
    "simulation.solve(t_eval, inputs={\"Current function [A]\": 1.6})\n",
    "\n",
    "# plot\n",
    "simulation.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PyBaMM can also simulate rest behaviour by setting the current function to zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd1dbdfee61f4e518c3534b44e0ce879",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=600.0, step=6.0), Output()), _dom_classes=('…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x79ae443f29f0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve the model at the given time points\n",
    "simulation.solve(t_eval, inputs={\"Current function [A]\": 0})\n",
    "\n",
    "# plot\n",
    "simulation.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading in current data <a name=\"data\"></a>\n",
    "\n",
    "To run drive cycles from data we can create an interpolant and pass it as the current function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Downloading file 'US06.csv' from 'https://github.com/pybamm-team/pybamm-data/releases/download/v1.0.0/US06.csv' to '/home/santa/.cache/pybamm'.\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd  # needed to read the csv data file\n",
    "\n",
    "model = pybamm.lithium_ion.DFN()\n",
    "\n",
    "# import drive cycle from file\n",
    "data_loader = pybamm.DataLoader()\n",
    "drive_cycle = pd.read_csv(\n",
    "    f\"{data_loader.get_data('US06.csv')}\", comment=\"#\", header=None\n",
    ").to_numpy()\n",
    "\n",
    "# load parameter values\n",
    "param = model.default_parameter_values\n",
    "\n",
    "# create interpolant - must be a function of *dimensional* time\n",
    "current_interpolant = pybamm.Interpolant(drive_cycle[:, 0], drive_cycle[:, 1], pybamm.t)\n",
    "\n",
    "# set drive cycle\n",
    "param[\"Current function [A]\"] = current_interpolant\n",
    "\n",
    "# set up simulation - for drive cycles we recommend using the CasadiSolver in \"fast\" mode\n",
    "simulation = pybamm.Simulation(model, parameter_values=param)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that when simulating drive cycles there is no need to pass a list of times at which to return the solution, the results are automatically returned at the time points in the data. If you would like the solution returned at times different to those in the data then you can pass an array of times `t_eval` to `solve` in the usual way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "513c59118a6144a8a205510aaa9e9527",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=600.0, step=6.0), Output()), _dom_classes=('…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x79ad7d74c620>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# simulate US06 drive cycle (duration 600 seconds)\n",
    "simulation.solve()\n",
    "\n",
    "# plot\n",
    "simulation.plot()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that some solvers try to evaluate the model equations at a very large value of `t` during the first step. This may raise a warning if the time requested by the solver is outside of the range of the data provided. However, this does not affect the solve since this large timestep is rejected by the solver, and a suitable shorter initial step is taken."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding your own current function <a name=\"function\"></a>\n",
    "\n",
    "A user defined current function can be passed to any model by specifying either a function or a set of data points for interpolation.\n",
    "\n",
    "For example, you may want to simulate a sinusoidal current with amplitude A and frequency omega. In order to do so you must first define the method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create user-defined function\n",
    "\n",
    "\n",
    "def my_fun(A, omega):\n",
    "    def current(t):\n",
    "        return A * pybamm.sin(2 * np.pi * omega * t)\n",
    "\n",
    "    return current"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the function returns a function which takes the input time.\n",
    "Then the model may be loaded and the \"Current function\" parameter updated to `my_fun` called with a specific value of `A` and `omega`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.lithium_ion.SPM()\n",
    "\n",
    "# load default parameter values\n",
    "param = model.default_parameter_values\n",
    "\n",
    "# set user defined current function\n",
    "A = model.param.I_typ\n",
    "omega = 0.1\n",
    "param[\"Current function [A]\"] = my_fun(A, omega)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that when `my_fun` is evaluated with `A` and `omega`, this creates a new function `current(t)` which can then be used in the expression tree. The model may then be solved in the usual way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "40a6ad712e184ba1a161932440d7a99d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=30.0, step=0.3), Output()), _dom_classes=('w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x79ad73fea6c0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set up simulation\n",
    "simulation = pybamm.Simulation(model, parameter_values=param)\n",
    "\n",
    "# Example: simulate for 30 seconds\n",
    "simulation_time = 30  # end time in seconds\n",
    "npts = int(50 * simulation_time * omega)  # need enough timesteps to resolve output\n",
    "t_eval = np.linspace(0, simulation_time, npts)\n",
    "solution = simulation.solve(t_eval)\n",
    "label = [f\"Frequency: {omega} Hz\"]\n",
    "\n",
    "# plot current and voltage\n",
    "output_variables = [\"Current [A]\", \"Voltage [V]\"]\n",
    "simulation.plot(output_variables, labels=label)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
      "[3] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[4] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.\n",
      "[5] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "612adcc456652826e82b485a1edaef831aa6d5abc680d008e93d513dd8724f14"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
